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Abstract 



We consider the application of finite-size scaling methods to isothermal- 
isobaric (constant-NpT) simulations of pure continuum fluids. A finite-size 
scaling ansatz is made for the dependence of the relevant scaling operators 
on the particle number. To test the proposed scaling form, constant pres- 
sure simulations of the Lennard-Jones fluid at its liquid-vapour critical point 
are performed. The critical scaling operator distributions are obtained and 
their scaling with particle number found to be consistent with the proposed 
behaviour. The forms of these scaling distributions are shown to be identi- 
cal to their Ising model counterparts. The relative merits of employing the 
constant-NpT and grand canonical (constant-// VT) ensembles for simulations 
of fluid criticality are also discussed. 
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I. INTRODUCTION 



The study of the phase behaviour of simple and complex fluids by computer simulation 
is a subject of much current research activity |]]. Of particular interest are the critical 
point properties of such systems @. In the vicinity of a critical point, the correlation 
length grows extremely large and may exceed the linear size of the simulated system. When 
this occurs, the singularities and discontinuities that characterise critical phenomena in the 
thermodynamic limit are shifted and smeared out P|. Unless care is exercised, such finite- 
size effects can lead to serious errors in computer simulation estimates of critical point 
parameters. 

To cope with these problems, finite-size scaling (FSS) techniques have been developed 
FSS methods enable one to extract accurate estimates of infinite- volume thermo- 
dynamic quantities from simulations of finite-sized systems. To date, their application to 
fluid criticality has been principally in conjunction with simulations in the constant-yuVT or 
grand-canonical ensemble (GCE). The principal merit of this ensemble is that the particle 
density fluctuates on the scale of the system as a whole, thus enabling direct measurement 
of the large-scale density fluctuations that are the essential feature of fluid criticality. The 
GCE has proven its worth in FSS studies of criticality in a variety of fluid systems including 
the Lennard- Jones (LJ) fluid [[||| and a 2D spin fluid model [J7]. 

Notwithstanding its wide utility, however, there exist many complex fluids for which use 
of the GCE ensemble is by no means efficient. Systems such as semi-dilute polymer solutions 
are difficult to simulate in the GCE due to excluded volume effects which hinder chain 
insertions. While smart insertion techniques go some way to ameliorating this difficulty, the 
long chain lengths of greatest interest are currently inaccessible 0. Similarly, electrolyte 
models such as the restricted primitive model show very poor acceptance rates for particle 
insertions due to the formation of bound ion clusters [[|. Thus it is interesting to ask whether 
one can deal with the near-critical density fluctuations in such systems without having to 
implement inefficient particle transfer operations. 

The approach we consider here, is to employ an ensemble wherein the total particle 
number is fixed, but the density is allowed to fluctuate by virtue of volume transitions. 
Specifically, we consider how the FSS ideas, hitherto only applied to systems with constant 
volume, may be generalised to an isothermal-isobaric (NpT-ensemble) simulation. Since 
finite-size scaling usually rests on the idea of comparing the correlation length with the 
(fixed) linear dimensions of the systems, the generalisation to systems whose linear dimen- 
sions are dynamically fluctuating is not completely obvious. We make a scaling ansatzior the 
near-critical scaling operator distributions and scaling fields, expressed in terms of powers of 
the particle number. This is then tested via a simulation study of the critical Lennard- Jones 
fluid, in which it found that the FSS predictions are indeed consistent with the simulation 
results. Finally we discuss the relative merits of the NpT- and /iVT- (GCE) ensembles for 
simulation studies of fluid criticality. 

II. FINITE-SIZE SCALING THEORY 

We consider a classical single component fluid, whose configurational energy (which we 
write in units of fcgT) resides in a sum of pairwise interactions amongst the N particles it 
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contains 

N 

*(M)= E 0(1^-^1) (i) 

i<j=i 

where </>(| — |) is the interparticle potential which, for this work, we assign the Lennard- 
Jones form: 

4>{r) = Aw[(a/r) 12 - (a/rf] (2) 

where w is the dimensionless well depth and a serves to set the length scale. 
Within the isothermal-isobaric ensemble, the partition function is given by 



N 



n{/m^ eh#(W)_pV1 ( 3 ) 



where p = P/ksT is the pressure, and V = L d is the homogeneously fluctuating system 
volume. The associated (Helmholtz) free energy is 

G = -k B T\nZ N = fc B T($({f}) +pV) (4) 

In the vicinity of the liquid- vapour critical point, the coarse-grained properties of the 
system are expected to exhibit scaling behaviour. For simple fluids with short-ranged in- 
teractions, this behaviour is Ising-like || and is controlled by two relevant scaling fields r 
and h that measure deviations from criticality. In general (in the absence of the special 
Ising 'particle-hole' symmetry), these scaling fields comprise linear combinations JlO] of the 
reduced chemical potential p and well depth w: 

t = w c — w + s(fi — fi c ), h = /i — /i c + r{w c — w) (5) 

where subscripts denote critical values, and s and r are non-universal field mixing parame- 
ters, the values of which depend on the specific system under consideration. 
The respective conjugate operators are defined by 

1 d\nZ N 1 d\nZ N 

{£) = ly)~dT~' {M} = (yj-efT (6) 

whereupon, utilizing equation |3|, and noting that dP = gdT + pdfi (with g = S/V, the 
entropy density), one finds 

1 1 

S = - [u-rp], M = - [p-su] (7) 

1 — sr 1 — sr 

where p = N/V is the particle density and u = &({r})/V is the energy density. We term 
Ai the ordering operator, while S is termed the energy-like operator. 

For the finite-size scaling behaviour of the joint distribution Pn(M.,£) we make the 
following ansatz: 

P N (M,S) ~ a^a^N-NypiaJ^N^M - M c ), a^ 1 ^^ - £ e ), a^iV 1 "^, aeN^r) (8) 
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where and ag are non-universal metric factors and P is a universal function in the sense 
that it is identical for all members of a given universality class and a given choice of boundary 
conditions. Here we have chosen simply the particle number N rather than the volume V as 
a measure of the 'finite-size' of the systems, using now suitable powers of N in the scaling 
assumption. 

The unspecified exponents x and y in equation |] can be related to the standard critical 
indices (3 and v by an argument analogous to that invoked by Binder |TT] for the Ising model 
in the canonical ensemble. Consider the generalised isothermal compressibility, k t , defined 
by 

This may be reexpressed as a fluctuation relation by appeal to equations |3] and |], whence 
one finds 

K T =f3(p)- l N((M-M c ) 2 ) (10) 

From equation [g, the scaling of the fluctuation in the ordering operator M. at finite r is 
given by 

{{M-M c ) 2 ) = aJ^N x J dM{M - M c ) 2 p N (aj}N x (M - M c ),N^r) (11) 

= N~ 2x a 2 M J dzz 2 p N (z,N 1 - y r) (12) 
= N~ 2x a 2 M /(iV 1 -^) (13) 

where f(N v r) is an unspecified scaling function. 

Now in the critical (r — > 0) and thermodynamic limits (N — > 00), one requires for 
consistency that k t ~ t~ 7 . Accordingly f(N y r) ~ (N 1 ~ y r)~' y , and matching powers of iV 
in equation [l(| and [13] then yields 

-2x- 7 (l-y) = l, (14) 
a relation that is satisfied by the assignments 

* = i (15) 

v = i-l, (16) 

where we have employed the hyperscaling relation dv = 7 + 2(3 to eliminate 7 from equa- 
tion |TJ. 

Thus precisely at criticality, where t = h = 0, equation § can be written: 

P N (M, £) ~ a£aJ 1 NW*N&- 1)/ *P*(ajJiN f,/ ' b '(M - M c ), a^N^- 1)/du (£ - £ c )), (17) 

where P* is a universal function characterising the Ising fixed point. We note that this 
scaling behaviour is simply related to that for the /iVT-ensemble by the substitution 
L d — ► iV. In what follows we will test the validity of equation [L7| via simulations of the 3D 
Lennard- Jones fluid at its liquid- vapour critical point. 
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III. MONTE-CARLO STUDIES 



Monte-Carlo simulations of the Lennard- Jones fluid were carried out within the NpT- 



ensemble using a Metropolis algorithm ||12|| . In common with many other simulation studies 
of the model, the potential was cutoff at a distance r c = 2.5<r, and was not shifted. In 
the course of the simulations, three separate periodic systems comprising N = 216, 363 and 
512 particles were studied. In order to minimise fluctuations in the acceptance rate for 
volume moves, the method of Eppenga and Frenkel [II] was employed. Within this scheme 
the system performs a random walk in In V rather than in V itself, thereby automatically 
adjusting the step size in the volume to be larger at lower densities than higher densities. 
This technique is particularly useful in the coexistence region region, due to the large volume 
fluctuations. 

From a previous detailed simulation study of the LJ fluid in the grand-canonical ensemble 
||, the critical temperature is known to lie at T c = 4/w c = 1.1876(3). Using this value of 
w, a number of short runs were performed in which the pressure p was tuned in order 
to locate approximately the coexistence curve. Once a near-critical value of p had been 
obtained, a long run comprising 6 x 10 6 volume change attempts and an equal number of 
displacements attempts per particle were performed. During this run, the histogram of the 
joint distribution P N (p,u) was accumulated. 

In order to locate more precisely the critical pressure, the ordering operator distribution 
Pn(A4) was studied for each value of N. Precisely at criticality, and for the correct choice 
of the field mixing parameter s, P^{M) is expected to be a symmetric function by virtue 
of the fact that the field mixing transformation restores the Ising symmetry of the critical 
fluid . The value of the field mixing parameter s and the pressure p were therefore tuned 



simultaneously using the histogram reweighting scheme [15|, until a symmetric function was 



obtained. This occurred for values s = —0.10(1) and p = 0.1093(6). This value of s is 
in accord with that obtained previously in the GCE study of the same model ||. The 
resulting form of Pn{M) for each system size, expressed in terms of the scaling variable 
x = aj}N^ du (M - M c ), is plotted in figure § The non-universal scale factor a^. was cho- 
sen so that the data for the iV = 512 system has unit variance, while the exponent ratio (3/v 
was assigned the Ising estimate (3/v = 0.518 of reference fL7 |. Also shown for comparison 
(solid line) is the scaled critical magnetisation distribution of the 3D Ising model obtained 
in a separate study fl~6 |. Clearly the data for each system size collapse well on to one an- 
other and are thus consistent with the proposed scaling form. They also agree well with the 
universal form of the Ising magnetisation distribution. We attribute the small deviations to 
corrections-to-scaling associated with the rather small numbers of particles. These discrep- 
ancies are comparable in size and form with those observed in reference JJ. Unfortunately, 
owing to these corrections-to-scaling, as well as to the somewhat poor statistical quality (a 
point to which we return in section |TV|), it was not possible to obtain a reliable independent 
estimate of the exponent ratio (3 jv. 

A procedure similar to that outlined above was also carried out for the energy-like op- 
erator P/v(£). In this case, however, there is no symmetry requirement for the critical form 
of the function |J. However, we found that by assigning the field mixing parameter r, the 
value r = —1.02 found in the previous GCE study of the LJ fluid 0, the function matched 
closely the form of the critical Ising model energy distribution. Again the scaling of the 
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Pn{£) with N was found to be consistent with the scaling form eq. 



IV. DISCUSSIONS AND CONCLUSIONS 

The aim of this work was to develop a FSS formalism for simulations of near-critical 
fluids at constant pressure. A scaling form was proposed and tested for the Lennard- Jones 
fluid, whose critical point parameters are known to high precision. Good consistency with the 
scaling prediction was found. Additionally, it was shown that the forms of the critical scaling 
operator distributions are the same as those for the /zVT-ensemble (and hence the canonical 
ensemble of the critical Ising model), despite the very different nature of the simulation 
ensembles. As is well known, finite-size effects may differ in various ensembles of statistical 
mechanics: the microcanonical ensemble, canonical and grand canonical ensembles of fluids 
are known to have distinct finite-size properties, and are characterized by different scaling 
functions. Our interpretation of the equivalence of the constant-NpT and constant-/z^T 
ensemble scaling functions, is that in the thermodynamic limit, p = p(T,p,), and therefore 
it is immaterial whether one uses fiT or pT as the pair of given intrinsic thermodynamic 
variables. All that matters is that both ensembles have a single extensive variable (N or 
V). Different scaling properties are obtained when two extensive variables are used, such 
as in the constant-NVT (canonical) ensemble of fluids. It follows that the GCE-based FSS 
techniques developed in reference |§ for accurately locating fluid critical point parameters, 
carry over directly to the NpT-ensemble. It is hoped that use of this latter ensemble will 
prove beneficial in situations where the GCE is highly inefficient, such as for long-chain 
polymer or electrolyte models. We add the further comment that finite-size scaling with 
N rather than V is well known in mean-field spin systems, when every spin interacts with 



every other spin, and geometry and space have no meaning [|T^ . 

While use of the NpT-ensemble is likely to be much more efficient than the GCE in cases 
where particle insertions have a low acceptance probability, we believe that for simple fluids, 
use of the grand canonical ensemble is much to be preferred. The present study of the LJ 
fluid consumed circa 300 hours CPU time on a Cray-YMP, using a vectorized program. By 
contrast, the previous GCE study of the same model was performed with considerably less 
computational effort using workstations. Moreover, it was possible to study much larger 
systems (comprising up to ~ 1800 particles), with considerably greater statistical accuracy 
than attainable here. The reasons for this difference in efficiency seem to be manifold. 
Owing to the fluctuating volume, it is not possible to easily implement a cell structure 
for efficiently locating neighbouring particles within the cutoff range. A calculation of all 
particle separations prior to a trial volume transition involves an 0(N 2 ) operation, which 
can only be avoided if no potential cutoff is employed [PHI - By contrast, use of a cell structure 



in the GCE leads to a O(N) growth in computational complexity. A second disadvantage 
of the NpT-ensemble, is that both volume changes and particle moves must be performed. 
For each particle move, two partial energy calculations are required, while for a volume 
change (with finite interaction cutoff), a total energy calculation is required. By contrast, 
only particle insertions need be implemented in a GCE scheme ||, each of which requires 
only a single partial energy calculation. Additionally it seems likely that for a given (p) 
the random walk in V required to pass from one phase to the other at coexistence is longer 
than that in N for the GCE, even if one employs a dynamically adjusting volume step size 
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such as that of reference |14|| . Consequently the correlation time in the NpT-ensemble is 
considerably greater than that of the GCE. 

Finally we note that another recent study has also attempted to apply finite-size scaling 
methods to the LJ fluid in the NpT-ensemble Hl9| . In this study the authors focused on the 
distribution of the specific volume v = V/N and attempted to locate the critical point by 
requiring scale invariance in Pn(v), subject to the requirement that Pn(v) has two peaks of 
equal height at coexistence. However, as we have shown in this work, the volume distribution 
is not the scaling counterpart for fluids of the Ising magnetisation, rather it is the operator 
M.. Moreover, the limiting (large N) critical point form of Pn{v) does not have two peaks 
of equal heights, and is in fact highly asymmetic in form, as figure shows. We therefore 
ascribe the rather low accuracy in the estimates for T c and the difficulties in obtaining the 
surface tension exponent in reference |19[ to this incorrect choice of the scaling operator. 
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FIG. 1. The critical scaling operator distributions Pn(A4) for the three system sizes 
N = 216, 363 and N = 512, expressed in terms of the scaling variable x 
The solid line is the critical magnetisation distribution of the 3D Ising model p 
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FIG. 2. The specific volume distribution Pn(v) for the N = 512 system, measured at the 
assigned values of the critical parameters. 
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